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METHODS AND APPARATUS FOR ADAPTIVE SIGNAL PROCESSING 
INVOLVING A KARHUNEN-LOEVE BASIS 

5 Field of the Invention 

The present invention relates generally to signal processing techniques and devices, and more 
particularly to signal processing techniques and devices which involve the utilization of a Karhunen- 
Loeve (KL) basis. 

10 Background of the Invention 

Many areas of signal processing utilize so-called eigendecompositions of covariance 
matrices. Such decompositions generally provide a basis, commonly known as a Karhunen-Loeve 
(KL) basis, in which signal expansion coefficients are uncorrected. The KL basis or a signal 
represented in the KL basis can immediately reveal crucial information in many applications. For 

15 example, in an array signal processing application, the KL basis is used to estimate directions of 
arrival of plane waves. This application allows the received signal to be separated into a signal 
subspace and a noise subspace, and is therefore called a subspace method. The signal space often 
changes with time, necessitating the use of subspace tracking techniques. 

Transform coding is another area in which it is useful to track a KL basis. In transform 

20 coding, a KL basis gives the best representation for encoding a Gaussian vector with scalar 
quantization and scalar entropy coding. When a KL basis is not known in advance or varies slowly, 
adaptive estimation or tracking can enhance performance. 

Gradient methods for adaptive signal-subspace and noise-subspace estimation are described 
in J.-F. Yang et al., "Adaptive Eigensubspace Algorithms for Direction or Frequency Estimation and 

25 Tracking," IEEE Trans. Acoust. Speech Signal Proc, Vol. 36, No. 2, pp. 241-251, February 1988. 
However, these methods are complicated both numerically and theoretically by orthonormalization 
steps. Adaptation of a Givens parameterized transform, which eliminates the need for explicit 
orthonormalization, was suggested as a method for finding a Karhunen-Loeve transform in P.A. 
Regalia, "An Adaptive Unit Norm Filter With Applications to Signal Analysis and Karhunen-Loeve 

30 Transformations," IEEE Trans. Circuits and Systems, Vol. 37, No. 5, pp. 646-649, May 1990. 
However, this approach fails to provide suitable convergence results. Subsequent work either does 
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not address step size selection, e.g., P.A. Regalia et al., "Rational Subspace Estimation Using 
Adaptive Lossless Filters," IEEE Trans. Signal Proa, Vol. 40, No. 10, pp. 2392-2405, October 1992, 
or considers only the most rigorous form of convergence in which step sizes must shrink to zero, 
gradually but not too quickly, e.g., J.-P. Delmas, "Performances Analysis of Parameterized Adaptive 
5 Eigensubspace Algorithms," Proc. IEEE Int. Conf. Acoust., Speech and Signal Proa, Detroit, MI, 
pp. 2056-2059, May 1995, J.-P. Delmas, "Adaptive Harmonic Jammer Canceler," IEEE Trans. 
Signal Proa, Vol. 43, No. 10, pp. 2323-2331, October 1995, and J.-P. Delmas, "Performances 
Analysis of a Givens Parameterized Adaptive Eigenspace Algorithm," Signal Proa, Vol. 68, No. 1 , 
pp. 87-105, July 1998. The latter form of adaptation, with step sizes approaching zero, is generally 

10 not suitable for use in practical applications. 

For details on other conventional basis tracking techniques, see J.-F. Yang et al, "Adaptive 
High-Resolution Algorithms For Tracking Nonstationery Sources Without the Estimation of Source 
Number," IEEE Trans. Signal Proa, Vol. 42, pp. 563-571, March 1994, P.A. Regalia, "An Unbiased 
Equation Error Identifier and Reduced-Order Approximations," IEEE Trans. Signal Proa, Vol. 42, 

15 No. 6, pp. 1397-1412, June 1994, B. Champagne, "Adaptive Eigendecomposition of Data 
Covariance Matrices Based on First-Order Perturbations," IEEE Trans. Signal Proa, Vol. 42, No. 
10, pp. 2758-2770, October 1994, B. Yang, "Projection Approximation Subspace Tracking," IEEE 
Trans. Signal Proa, Vol. 43, No. 1, pp. 95-107, January 1995, and B. Champagne et al., "Plane 
Rotation-Based EVD Updating Schemes For Efficient Subspace Tracking," IEEE Trans. Signal 

20 Proa, Vol. 46, No. 7, pp. 1886-1900, July 1998. 

The above-identified conventional techniques for tracking a KL basis fail to provide adequate 
performance, such as local convergence within specified step size bounds. In addition to exhibiting 
a lack of suitable convergence guarantees, most conventional techniques are computationally 
complicated and not highly parallelizable. A need therefore exists for improved techniques for 

25 tracking a KL basis. 
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Summary of the Invention 

The present invention provides improved signal processing techniques and devices for 
Karhunen-Loeve (KL) basis tracking through the use of stochastic approximations to specified 
gradient descents. In order to keep computational requirements low and to provide a degree of 
5 parallelizability, the signal processing techniques and devices of the present invention utilize 
computations which have locality and regularity properties. 

In accordance with the invention, a transform is represented in a reduced-parameter form, 
such as a Givens parameterized form or a Householder form, such that the reduced-parameter form 
for mNx N transform comprises fewer than TV 2 parameters. An updating process for the transform 
10 is implemented using computations involving the reduced-parameter form, and an adaptation of the 
transform is represented directly as one or more changes in the reduced-parameter form. 

An illustrative embodiment of the invention provides a gradient descent algorithm 
particularly well suited for transform adaptation in transform coding applications. The algorithm 
can be implemented in one or more signal processing devices, such as a digital signal processor 
15 (DSP), filter, encoder, decoder, etc. The algorithm in the illustrative embodiment uses a pairwise 
energy compaction property of the KL transform. Advantageously, the algorithm converges locally 
in mean within specified step size bounds. In addition, the algorithm allows trade-offs to be made 
between speed of convergence and steady-state error. It is particularly simple to implement and 
exhibits a rate of convergence which is generally better than that of more computationally 
20 demanding conventional approaches. 

The invention provides a number of other advantages over conventional techniques. For 
example, the invention avoids the need for an explicit eigendecomposition operation in 
implementing the transform, and does not require explicit estimation of the autocorrelation matrix 
of the source signal being processed. In addition, the invention provides better parallelizability of 
25 computations than the previously-noted conventional techniques, and allows the transform update 
process to be driven by quantized data. 

The algorithm in the illustrative embodiment preferably operates on Givens parameterized 
matrices, i.e., matrices which are a product of multiple Givens rotations, so as to ensure 
orthogonality, increase numerical robustness, and minimize the number of parameters. 
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Although particularly well suited for use in transform coding applications, the invention can 
also be applied to other types of signal processing applications, including array signal processing, 
general-purpose signal filtering, or other applications involving a KL basis. Moreover, the 
techniques of the invention are applicable to a wide variety of different types of signals, including 
5 data, speech, audio, images, video and other types of signals. 



Brief Description of the Drawings 

FIG. 1 is a schematic representation of a signal processing device which implements a Givens 
parameterized form of an orthogonal transform in accordance with the invention. 
10 FIG. 2 is a schematic representation of a signal processing device which implements a 

stochastic approximation of a gradient descent algorithm in a Givens parameterized form in 
accordance with the invention. 

FIGS. 3 A and 3B show plots of deterministic simulations which demonstrate differences in 
rates of mean convergence for gradient descent algorithms as applied to example source signals. 
15 FIGS. 4 A, 4B arid 4C show plots of stochastic simulations for gradient descent algorithms. 

FIG. 5 shows a schematic representation of a signal processing device for backward adaptive 
transform coding in accordance with the invention. 

FIG. 6 shows plots of simulations of a backward adaptive version of a gradient descent 
algorithm in accordance with the invention. 

20 

Detailed Description of the Invention 

The invention will be illustrated below in conjunction with exemplary signal processing 
techniques and devices. The techniques and devices described may be applied to processing of a 
wide variety of different types of signals, including data signals, speech signals, audio signals, image 
25 signals, and video signals, in either compressed or uncompressed formats. The term "vector" as used 
herein is intended to include any grouping of coefficients or other components representative of at 
least a portion of a signal. The term "signal processing device" is intended to include any type of 
device suitable for providing at least a subset of the processing operations of the present invention, 
such as a digital signal processor (DSP), microprocessor, computer, application-specific integrated 

4 



V.K. Goyal 5 

circuit (ASIC), filter, encoder, decoder, etc. The particular device implementation used in a given 
application will of course be a function of the particular requirements of that application. 

As described previously, it is often desirable in signal processing applications to update a 
Karhunen-Loeve (KL) basis estimate as data vectors are processed. An illustrative embodiment of 
the invention to be described below provides a Givens parameterized gradient descent algorithm for 
updating a KL basis estimate. This embodiment, which is based on pairwise maximization of coding 
gain, is contrasted herein to two conventional algorithms based on stochastic gradient descents of 
global cost functions. Although all three of these algorithms can be shown to converge locally in 
mean, the algorithm of the illustrative embodiment of the invention converges much faster than the 
other two while also being computationally the simplest of the three. It should be emphasized, 
however, that the illustrative algorithm is presented by way of example, and should not be construed 
as limiting the scope of the invention. 

The three algorithms to be described herein are also referred to as Algorithm 1, Algorithm 
2 and Algorithm 3. It should be understood that Algorithms 1 and 2 are the conventional algorithms 
based on stochastic gradient descents of global cost functions, and are described for purposes of 
comparison with Algorithm 3, which is the illustrative embodiment of the invention based on 
pairwise maximization of coding gain. 

The approach utilized in Algorithms 1 , 2 and 3 may be analogized in some respects to Wiener 
filtering, as described in, e.g., P.M. Clarkson, "Optimal and Adaptive Signal Processing," CRC 
Press, Boca Raton, FL, 1993, which is incorporated by reference herein. For optimal 
finite-impulse-response Wiener filtering, one requires knowledge of second-order moments of 
signals and one needs to solve a linear system of equations, which in some contexts may be 
considered expensive. The linear system of equations is commonly replaced with a gradient descent 
and the exact moments are replaced with immediate stochastic estimates, leading to the well-known 
least-mean-squares (LMS) algorithm. By analogy, determining the optimal transform for transform 
coding requires the determination of an eigendecomposition, typically an expensive computation. 
Algorithms 1 , 2 and 3 replace this computation with a gradient descent. It should also be noted that 
both the LMS algorithm and the gradient descent algorithms can be driven by quantized data. 
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The following notation will be used in describing Algorithms 1, 2 and 3. Let {x„} neN be a 
sequence of R N - valued random vectors and let X n = E [ x n x T n ]. It is assumed for simplicity and 

clarity of illustration that the dependence ofX n on n is mild. The algorithms are designed to produce 
a sequence of orthogonal transforms T n such that T n X n T n is approximately diagonal for each n. The 

5 algorithms are preferably strictly causal, i.e., T n should depend only on \x k j k o . In general, X n is 

not known, but instead must be inferred or otherwise estimated from \x k \ k _ Q • 

An orthogonal T n such that T n X n T T n is diagonal is called a Karhunen-Loeve transform (KLT) 

of the source vector jt n . Finding a KLT is equivalent to finding an orthonormal set of eigenvectors 
of the symmetric, positive semidefinite covariance matrix X n ; i.e., solving a so-called symmetric 

10 eigenproblem, as described in G.H. Golub et al., "Matrix Computations," Johns Hopkins Univ. 
Press, Baltimore, MD, Second Edition, 1989, which is incorporated by reference herein. Thus, 
finding an algorithm for transform adaptation may be viewed as an attempt to approximately solve 
a sequence of symmetric eigenvalue problems with noisy data. 

The present invention avoids the need for an explicit eigendecomposition by using incoming 

1 5 data samples to make small adjustments in the transform so that changes in the KL basis are tracked 
and any initial misadjustment is eliminated. Algorithms 1 and 2 make the transform adjustments in 
the direction of the negative gradient of a cost function that is minimized by the desired 
diagonalizing transform. Examples of such cost functions and corresponding gradient descent 
analysis with respect to these functions will be given below. Algorithm 3 applies a gradient descent 

20 to an energy compaction property that holds for any pair of transform coefficients, instead of 
minimizing a single global cost function through gradient descent. 

Certain additional notation used in the description of Algorithms 1, 2 and 3 will now be 
introduced. A Givens rotation G of # radians counterclockwise in the (i 9 j) coordinate plane is an 
identity matrix except for 

25 
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This rotation will be denoted G ij6 . A Givens rotation is also sometimes referred to as a Jacobi 
rotation. 

5 Any N x N orthogonal matrix can be factored as a product of AT = N(N- 1 )/2 Givens rotations 

with i < j and each rotation angle 0e [-7T/4, tu/4). See, for example, F. J. Vanpoucke et al., "Factored 
Orthogonal Transformations For Recursive Eigendecomposition," IEEE Trans. Circuits Syst.-II: 
Anal. Dig. Sig. Proa, Vol. 44, No. 3, pp. 253-256, March 1997, which is incorporated by reference 
herein. 

1 0 Givens parameterizations of orthogonal matrices will be used extensively in the description, 

so the following notation that eliminates one subscript will be used. To index only the needed (/, 
j) pairs, map (z,y) to where (/,;) is the Ath entry in a lexicographical list of (ij) e{l,2, . . . , N} 2 
pairs with i < j. The reverse of this index mapping is denoted by (i k J k ) being the pair corresponding 
to index k. A Givens parameterized representation is given by 

15 



where 0 = (0 l9 6 2 , . . . , 0 K ) and the indices have been remapped according to G k B = G ikJk $ . 

FIG. 1 shows a signal processing device which implements the above-described Givens 
20 parameterized form of an orthogonal transform with the index remapping for a case of 7V= 4. Each 
G k block in the figure represents a Givens rotation matrix. The Givens parameterized form ensures 
orthogonality and minimizes the number of parameters. Matrices in this form can be multiplied 
using conventional techniques such as those described in, e.g., the above-cited FJ. Vanpoucke et al. 
reference. This is an important aspect of the present invention, as it produces a multiplicative update 
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matrix in Givens parameterized form, i.e., an orthogonal matrix U Qn in Givens form is generated 
such that T n+l = U Qn -T n 

Global Gradient Descent Algorithms (Algorithms 1 and 2) 

5 This section will consider diagonalizing a fixed matrix X 9 so its time index is not shown. The 

problem of finding a KLT T corresponding to a co variance matrix Xcan be cast as a minimization 
problem if there is a cost function that is minimized if and only if T is a KLT. Many such cost 
functions exist, and two are described in greater detail below, along with corresponding analyses of 
deterministic gradient descents of these cost functions. These analyses lead directly to 
10 convergence-in-mean results for the stochastic transform update algorithms, i.e., Algorithms 1 and 
2. 

One type of example cost function is the squared 2-norm of the off-diagonal elements of Y 

= TXT T : 

N N 
;=1 j=\J*i 

This cost function attains its minimum value of zero if and only if T exactly diagonalizes X. 
Another example of a cost function is 

20 ^ 2 01=fk, 

l=\ 

which, although less clearly connected to the diagonalization of X, is also minimum when T 
diagonalizes X. This can be seen as follows. As T varies over the set of orthogonal matrices, the 
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diagonal elements of Y lie in the convex hull of the permutations of the eigenvalues of X, as 
described in A.W. Marshall et al., "Inequalities: Theory of Majorizations and Its Applications," 
Academic Press, San Diego, CA, 1979, which is incorporated by reference herein. Furthermore, the 
sum of the diagonal elements of Y is constant. The minimum of the product occurs when the 
elements are as far apart as possible, i.e., at a corner of the convex set. Since the corners are 
permutations of the eigenvalues, J 2 is minimized when Tis a diagonalizing transform. 

In high resolution transform coding with optimal bit allocation, the mean-squared error 



(MSE) distortion per component at rate R bits per component is proportional to p 2 (Y)2- 2R . In 

this context, it is shown in A. Gersho et al., "Vector Quantization and Signal Compression," Kluwer 
1 0 Acad. Pub., Boston, MA, 1 992, which is incorporated by reference herein, that J 2 is minimized when 
TisaKLT. 

Both of the above-described cost functions J x andJ 2 are continuous in each component of Y. 
It should be noted that many other qualitatively different cost functions may be used, including, e.g., 
non-negative-coefficient linear combinations of the cost functions J x and J v 
15 Given a cost function J(Y), defined on symmetric positive semidefinite matrices, that is 

minimized if and only if Y is diagonal, one skilled in the art of adaptive filtering could approach the 
design of an adaptive KL transform as follows. First, compute the gradient 



V T J(TXT T ) , 



and then update according to 



T n+l = T n -»V T J(TXT T ) 



T=T„ 



25 where u is a small, positive step size. The initial difficulty with such an approach is that the 
orthogonality of T n does not imply the orthogonality of T n+1 . Furthermore, the gradient will generally 
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be dense, making computations difficult. Using a descent over the Givens parameterized form of 
the transform alleviates the first difficulty, but not the second. 

For clarity of illustration, deterministic descents using the cost functions J x and J 2 will now 
be described. Using the Givens parameterized representation of the transform matrix, finding a 

5 diagonalizing transform amounts to minimizing J x or J 2 over 0 G [-71/4,71 /A)" . A straightforward 
gradient descent would update the parameter vector through 

© w+ i = 0„-^v 0 4c/ 0 ^/ 0 r )| 0=0n5 (i) 

10 with step size // e M + . This approach can work, but is not computationally attractive because the 
gradient is dense, i.e., each component of V J i depends on every component of 0. As an alternative, 

T 

gradient descent can be used to compute a multiplicative update matrix. With Y n = T n XT n , one can 
update by 



15 T n+l = U @ T n , (2) 



where 



0„ = -^V 9 J,(c/ 0 rX)|e=o (3) 

20 

with step size u e R + . As will now be detailed below, the gradient in equation (3) is sparse and easy 
to compute for both cost functions J X (Y) and J 2 (Y). 

For notational convenience, define A (k) , 1 < k < K, elementwise by: 
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4 k) = dZ v /dB k .. 



where 



Z=U e YU T e 



Expanding C/^and using the product rule of differentiation gives 



A w = GA • • • G k _,H k G M ■ ■ ■ G K YG T K ■ ■ ■ G T 2 G T X + 



oik) 



G { G 2 • • • G K YG T K • • G T M H T k G T k _ x • • • G 2 G[ 



o(k)T 



(4) 



10 where H k is the derivative with respect to 0 k of the Givens matrix i.e., H k is all zeros except for 



- sin9 k cosG k 
-cosG k - sin9 k 



Evaluating at 0 = 0 makes all the G, identity matrices. H k is sparse, giving 



15 



B 



(k) 



0 = 0 



Y J,2 



7* 



(5) 
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where only the nonzero entries, in the i k and j k rows, are shown. 

As another notational aid, define ^ ( *'° elementwise by Af j) = Calculation of 

A {k,!) is simple and similar to that of A {k) and thus is not further described herein. For additional 
details, see V.K. Goyal, "Beyond Traditional Transform Coding," Ph.D. Thesis, Univ. California, 
Berkeley, 1998, Published as Univ. California, Berkeley, Electron. Res. Lab. Memo. No. UCB/ERL 
M99/2, January 1999, which is incorporated by reference herein. 

The minimization of J x by deterministic descent will now be described. Using equation (5), 
the components of the gradient of J x can be expressed as 



dJ x (u e YUl) „ „ / (i) 

pi A ^ij^ij Lu ^tj\ ij ji j 



10 " B °' = 2Z l Z v A*> = 2Z t Z ll [ltf> + BK>)- (6) 

k i*j i*j 



The partial derivative (6) generally depends on all of the components of 0, but evaluating at 0 = 0 
gives a very simple expression. Cancellation of like terms yields 

= 4Y (y - Y ) n\ 

'kJk\ JkJk 'k'k J • w 

0 = 0 



15 



djSu B Yuj) 



Thus the gradient computation in equation (3) requires a constant number of operations per 
parameter, independent of N. It should be noted that since the Frobenius norm is unchanged by 

orthogonal similarity transformations, minimizing J x is equivalent to maximizing J X (Y)= ^. =1 Y* . 

20 Of course, J { + J\ is a constant, so V J x = - V J x , demonstrating an easy way to obtain equation 
(7) and that a gradient ascent of J[ is equivalent to the above-described gradient descent. 
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The equilibrium points of a descent with respect to J x are those for which equation (7) is zero 
for all k. The desired, stable equilibria are obtained when all Y t j are zero. The following gives step 

size bounds for local convergence around these equilibria. 

Denote the eigenvalues of Xby A, > A 2 > — A N . It can be shown that the gradient descent 
5 update algorithm described by equations (2), (3), and (7) converges locally to a diagonalizing 
transform if 



0< n < 



2max(l] - A, y ) 



2(X, - X N ) 



10 This convergence property of the illustrative gradient descent algorithm can be proven as 

follows. Because the interest is in local properties, it can be assumed without loss of generality that 
7= diag ([X x X 2 — A N ]). The key to the proof is linearizing the autonomous, nonlinear, discrete-time 
dynamical system described by equations (2) and (3). More particularly, describe the system as 

15 ' ®„ + i=0„-Meu 

which upon linearization about the desired equilibrium 0 gives 



0„ +1 = (/-^)0„, where F y = 




20 Since / is continuously differentiable, a sufficient condition for local convergence is that the 
eigenvalues of I - fiF lie inside the unit circle. See M. Vidyasagar, "Nonlinear Systems Analysis," 
Prentice-Hall, Englewood Cliffs, NJ, 1978, which is incorporated by reference herein. Evaluation 
of F then proceeds as follows. Differentiating equation (6) gives 



13 




d 2 J,(u 9 YUl) 



= 2Z 



(z s ^') + 4« 4 »). 



V.K. Goyal 5 



(8) 



F is obtained by evaluating (8) at 0 = 0. Z reduces to Y, which is diagonal, so the first term makes 
no contribution. Simplifying the second term using equations (4) and (5) gives 



59,56; 



©=o 



4(1, -i J nk-i 

0 otherwise 



(9) 



The eigenvalues of / - juF are 1 - 4\i[l ik -^ Jk ) - The proof is completed by requiring that these 

eigenvalues all lie inside the unit circle. 
1 0 The minimization of J 2 by deterministic descent will now be described. The gradient of J 2 

is given elementwise by 

djAujTUl) ( T \^( I (k) ) 



15 To evaluate at 0 = 0, first note that 



A {k) \ = 2B {k) \ 

mm I 0=0 mm \q-q 



2Y hm Xm=i k , 
0 otherwise. 
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Thus only the m = i k and m = j k terms of the sum contribute, and 



dJ 2 {u Q YUl) 




hJk \ JtJt hh 



(12) 



0=0 



5 The derivative (12) differs from that of (7) only by a factor of 2* 1 Y~J k YJ k ) k J~[ Y ii . Thus, the 

equilibria are the same except that any Y H = 0 gives an equilibrium. The occurrence of Y ti = 0 is a 
degenerate condition that would be accompanied by a row and column of zeros. It can only occur 
if Xhas a zero eigenvalue, meaning that the source vectors lie in a proper subspace with probability 
one. As in the case of minimization of J x by deterministic gradient descent, one can also prove a 
10 local convergence result for the desired equilibria in the above-described minimization of J 2 , as 
follows. 

Denote the eigenvalues of Xhy X x > X 2 > — > \ N . The gradient descent described by 
equations (2), (3), and (12) converges locally to a diagonalizing transform if 



This convergence result can again be shown by linearizing the autonomous, nonlinear, 
discrete-time dynamical system implicitly defined. The result is summarized as 



15 



0< \x < 






if k = /, 
otherwise. 
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Requiring the eigenvalues of / - \iF to lie inside the unit circle completes the proof. 

Stochastic descents will now be described for the two example cost functions J x and J 2 . For 
the Karhunen-Loeve filtering of a random signal, it is generally desirable to provide stochastic 
versions of the deterministic gradient descent algorithms previously described. More particularly, 
5 one would like to compute an update vector based on a single data vector such that the expected 
value of the update is the same as the update computed in the deterministic algorithm. This gives 
a stochastic algorithm that converges in mean under the step size conditions of the corresponding 
deterministic algorithm. 

In well-known heuristic analyses of the LMS algorithm, e.g., from the above-cited reference 
10 P.M. Clarkson, "Optimal and Adaptive Signal Processing," CRC Press, Boca Raton, FL, 1993, an 
"independence assumption" is used to simplify computations. Analyses based on the independence 
assumption are qualitatively accurate. One of the analyses below uses a similar independence 
assumption. 

The minimization of J x by stochastic gradient descent will now be described. This process 
1 5 is referred to herein as Algorithm 1 . To obtain an update following equation (7), let 

where the time index n has been suppressed. For Gaussian data, 

20 

Therefore, the expected parameter vector changes in an algorithm using equations (2) and (13) are 
the same as those in the deterministic iteration of equations (2) and (3) with (7). The factor of 3 
25 difference between equations (13) and (14) is specific to a Gaussian source. For a non-Gaussian 
source, one can absorb this factor into the step size selection. 
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Consider the adaptive filtering of a Gaussian signal with covariance matrix Xusing equations 
(2) and (13). It can be shown that the mean of the transform, in Givens parameterized form, 
converges locally under the conditions given in the above description of the corresponding 
deterministic descent. 

5 Notice that in equation (13) the 0 k parameter depends only on the i k and j k components of y. 

This is desirable because it allows a parallelization of computations and a regularity of data flow, 
as is illustrated below in conjunction with FIG. 2. 

The minimization of J 2 by stochastic gradient descent will now be described. This process 
is also referred to herein as Algorithm 2. In order to obtain such an algorithm, following equation 
10 (12), one possible choice is 



= -n 



n 



yi 



(15) 



If the source is Gaussian and, for the purpose of examining the kth parameter, it is assumed 
15 that y t and y h are independent of the remaining components ofy, 



2Y (y -y) 

*kJk \ JkJk ) 



(16) 



Consider the adaptive filtering of a Gaussian signal with covariance matrixXusing equations 
20 (2) and ( 1 5). Using the independence assumption that yields equation (16), it can be shown that the 
mean of the transform, in Givens parameterized form, converges locally under the conditions given 
previously in conjunction with the corresponding deterministic descent. 
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With the exception of the Y\ " ^ y] factor, which may be incorporated in an update calculation 

block 100 such as that shown in FIG. 2, the update rule (15) shares the desirable regularity and 
parallelizability features of the previously-described Algorithm 1 . 

FIG. 2 shows a schematic representation of signal processing device which implements a 
5 stochastic gradient update in a Givens parameterized form. The "V" blocks represent gradient 
computations based on equations (13) or (15). It should be noted that all data dependencies are not 
shown for equation (15) in this particular representation. The update calculation block 100 is a 
scaling of the gradient by step size ju, and the multiplication of Givens parameterized matrices may 
be implemented using conventional techniques, such as those described in, e.g., F.J. Vanpoucke, et 
10 al, "Factored Orthogonal Transformations For Recursive Eigendecomposition," IEEE Trans. 
Circuits Syst.-II: Anal. Dig. Sig. Proc, Vol. 44, No. 3, pp. 253-256, March 1997, which is 
incorporated by reference herein. 

Pairwise Energy Compaction Algorithm (Algorithm 3) 

1 5 Algorithms 1 and 2 in the previous section are based on a global property. More specifically, 

they are based on the diagonalizing property of a KLT; i.e., that the output of a KLT has a diagonal 
covariance matrix. In this section a different property of the KLT, i.e., maximal energy compaction, 
is used to derive a gradient descent algorithm illustrative of the present invention. 

With the conventional sorting of basis vectors, the output of a KLT satisfies the following 

20 property of maximizing the energy in the lowest-numbered components: 

m m 

Yy,,= max Yj Y u> for m=\,2,...,N. 

orthogonal T 

See, e.g., K.R. Rao, "Discrete Cosine Transform: Algorithms, Advantages, Applications," Academic 
25 Press, 1990, which is incorporated by reference herein. This property will be utilized in the 
following equivalent, but less common form 
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£lj= min for m=l,2,...,tf. (17) 

/=m orthogonal T /=m 



10 



15 



This form uses gradient descent instead of the gradient ascent of the previous form. Because of the 
essential uniqueness of the KLT, the energy packing property can be indirectly captured by the 
global cost functions 7, and J 2 . However, this section will not be concerned with global behavior, 
but will instead determine how the energy packing property manifests itself between pairs of 
transform coefficients. Having local optimality in each pair gives the desired global property as a 
byproduct. 

Referring again to the schematic representation of FIG. 1, if the transform is a KLT, then^ 4 
has minimum possible power, y 3 has minimum power given that>> 4 is fixed, etc. This suggests that 
a KLT can be obtained by adjusting each Givens rotation in FIG. 1 to minimize the power sent to 
the lower branch. As will be demonstrated below, this approach in fact leads to a locally convergent 
adaptation algorithm. 

Regardless of the size of the overall transform, the operation of rotation block G k is described 

by 











x" 









cos 8 k sin 0 k 
-sin8 k cos 0 k 



x' ik COsQ k + x^sin0 k 

-x] sin0 . + x' cos0 , 

i k Jk K 
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where jc' and x" denote the input and output, respectively, of that transform stage. Now since 



d(x'') 2 dx'' i x 

\ A * ; = 2x" = 2x7 \-x' cos0 k - x) sinG k = -2x7 x\' , (18) 
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increasing Q k by 2px" x" is an instantaneous gradient descent of the energy in the lower branch. 

This gradient is very convenient because it is simply the product of quantities that are computed in 
the transform itself. 

To ascertain the convergence in mean of the new iteration described above, the expected 
update 0£, derived from equation (18), is used: 



2n 



= £[(x£cos0 t + x' Jt sine sinG, + x' Jk cosQ k ) 



2 V 



jk 



- E 



[*; 2 ])sin2e t + 2< 



x\ x' 

h Jk 



cos20 



(19) 



10 Since x\ and x' Jt depend on all earlier rotation stages, the components of 0 are coupled in a 

complicated way. However, using the approach of the previous section, the descent can be used to 
compute a multiplicative update, as illustrated in FIG. 2. The resulting update algorithm is referred 
to herein as Algorithm 3. The convergence in mean then depends only on the behavior of the 
expected gradient (19) around 0. 
15 Denote the eigenvalues of Xby X } > X 2 > - > X N . It can be shown that the deterministic 

iteration described by equation (2) with 



hJk 



20 converges locally to a diagonalizing transform if 
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0< u < 



max(x. ( - A, j.j 



\ x -'k N 



Under the same step size bounds, it can be shown that the stochastic iteration described by equation 
(2) with 



0* = viy ik y 



(20) 



10 



gives a transform sequence that, in Givens parameterized form, converges locally. 

As in the proofs of convergence given for the deterministic descents described previously, 
the expected value of the parameter vector is described by a discrete-time dynamical system. 
Linearizing about the desired equilibrium, the stability of 



Q n+i = (l+ M F)® n where F tJ = 



d 

d6 s 



0=0 



is established. This yields for the diagonal elements of F, 
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d 


[ 9 *1 








] 


_ae/ 


. i 1 . 





= [2(4tf]- 4< 2 ])cos20 k - 44,>; t ]sin20,]| 0=o 

0-0 



(21) 



and the off-diagonal elements are 
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F kl = 



50, 



0=0 



5 (4<H< 2 ]) 

50, 



sm2Q k +2- 

->Oas0->O 



dE\x\ x' 

\_ l k Jk 



59, 



cos20 t 



= 0 . (22) 



»Oas0-»O 



0=0 



The eigenvalues of F are2|A,^ -X.J, £ = \,2,...,K; thus the eigenvalues of / + juF 

are 1 + // l{x jk - \ ) • When all the eigenvalues lie in the unit circle, the dynamical system converges 
5 locally to the desired equilibrium. 

Simulations and Comparisons 

The mean convergence rates of the above-described illustrative update algorithms will now 
be examined in greater detail. The linearizations used in the above convergence proofs facilitate 
10 estimation of rates of convergence. Consider Algorithm 1, the stochastic gradient descent with 
respect to J v In the vicinity of the desired equilibrium, the error in 6 k at the nth iteration is given 
approximately by 
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based on equation (9). For local stability in this example, the bracketed quantity is required to be 
in the interval (-1,1) for all k. Furthermore, the bracketed quantity should be as close to zero as 
possible, in order to provide fast convergence. The ability to do this through the choice of fi is 
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limited by the variability of (x k - X h j . Accordingly, one can define pseudo-eigenvalue spreads 



for the three illustrative algorithms as follows: 



max(\ t - X ,) 

Xx) ^ 



min(x i - X ) 

i<j \ J / 



2 5 ^2 



(JT) = " , *'\ 2 , and s 3 U)= ^ 



min- 

'</ AyA ■ 



max! A, . - A, ■ ) 
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These pseudo-eigenvalue spreads may be viewed as analogous to the eigenvalue spread in LMS 
adaptive filtering as described in the above-cited P.M. Clarkson reference, "Optimal and Adaptive 
Signal Processing," CRC Press, Boca Raton, FL, 1993. 

Notice that a pseudo-eigenvalue spread is not less than one and that s x (X) = (s 3 (X)) 2 zs 3 (X). 
This suggests that with appropriate step size choices, Algorithm 3 will have faster local convergence 
than Algorithm 1 for any signal The difference between s x (X) and s 2 (X) suggests that the superior 
choice between Algorithms 1 and 2 will depend on along with the choice of ju . These observations 
are confirmed through the following example. 

Consider the matrices 



X x = diag 



7 5 1 

l5 8'8'2 



and X 2 = diag 



1, 



III 
2'4'8 



15 for which s^X^ 16< 28= s 2 (x,)and s Y (x 2 ) = 49> = s 2 (x 2 ). Based on these pseudo-eigenvalue 



spreads, it is expected that the Algorithm 1 global descent with respect to J x will perform better than 
the Algorithm 2 global descent with respect to J 2 for matrix X l9 and vice versa for matrix X v For 



23 




V.K. Goyal 5 

Algorithm 3, s 3 (X x ) = 4 and s 3 (X 2 ) = 7. These are the smallest spreads, so it is expected that 
Algorithm 3 will have the fastest convergence for both of the matrices in this example. 

FIGS. 3 A and 3B show deterministic simulation results for each of the three algorithms for 
the respective X x andX, matrices of the above example. Algorithm 1 is the global gradient descent 
5 with respect to J x as in equation (13), Algorithm 2 is the global descent with respect to J 2 as in 
equation (15), and Algorithm 3 is the pairwise energy compaction descent in accordance with 
equation (20). The simulation results indicate that the performance is as predicted. In the 
simulations, each algorithm was used with the step size |i that maximized the speed of convergence 
of the slowest mode, as given in Table 1 below. Since the signal statistics are not known in advance, 
10 the step size selection will generally not be optimal. However, the pseudo-eigenvalue spread 
provides a bound on performance. 







X 2 


Algorithm 1 


32/17 


16/25 


Algorithm 2 


1024/145 


512/53 


Algorithm 3 


8/5 


1 



Table 1: Step Sizes Used in Simulations of Mean Convergence 



20 The simulation results are averaged over 1 00 randomly chosen initial transforms. Each curve 

in FIGS. 3 A and 3B thus represents the average of 100 simulations with randomly chosen initial 
transforms corresponding to 0 uniformly distributed on [-tt/4, ti/4). Convergence could also be 
measured by J, or J 2 , and results of this type are given in the above-cited Ph.D. thesis of V.K. Goyal. 
At each iteration, the error in the transform is measured by the Frobenius norm of the difference 

25 between a true KLT, i.e., an identity matrix, and the adapted transform. Since the interest is in 
finding any KLT, this approach uses the minimum norm over all row permutations and 
multiplications of rows by ±1 : 
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mm 

permutations P 
diagonal a s.t.o 2 = I 



I-aPT n 



Stochastic simulations have also been performed on the three illustrative update algorithms. 
In the stochastic setting, e.g., using equations (13), (15) or (20), transform updates are on average 
in the correct direction, but may have high variance. The step size should not actually be chosen 
5 close to the previously-described predicted bound for stability. The simulations described below in 
conjunction with FIGS. 4 A, 4B and 4C indicate that a good choice is one or two orders of magnitude 
less than the bound. 

Consider by way of example a zero-mean Gaussian source with covariance matrix 



1 0.6 0.36 0.218 

0.6 1 0.6 0.36 

0.36 0.6 1 0.6 

0.218 0.36 0.6 1 



FIGS. 4 A, 4B and 4C show stochastic simulation results for Algorithm 1, Algorithm 2 and 
Algorithm 3, respectively, using the example four-dimensional source X given above. The source 
is Gaussian with E[xpc] = 0.6 |,vl . In all simulations the initial transform was the identity transform. 
After finding the eigenvalues of X, predicted bounds for stability 

15 \x] « 0.1131, \i* 2 * 0.5864, and (i 3 « 0.4757 were computed for the three algorithms, respectively. 

The simulations used \i = |i7y for various values of y 9 where n* is the maximum step size for 
convergence in mean. The resulting curves are labeled by y in the figures. Each curve represents 
the average of 100 simulations. 

Based on intuition from LMS filtering, one would expect a trade-off between speed of 
20 convergence and steady-state error. The results shown in FIGS. 4 A, 4B and 4C confirm this 
intuition. More particularly, it can be seen from the curves that with a large step size (small y), 
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steady state is approached quickly, but the steady-state error is large, while with smaller step sizes 
(larger y), steady-state error is smaller but convergence is slower. In each of FIGS. 4A, 4B and 4C, 
the lower hull of the curves gives the minimum error for a given number of iterations, assuming a 
fixed step size. This lower hull is lowest for Algorithm 3. The superiority of this algorithm is also 
5 predicted by comparing pseudo-eigenvalue spreads. 

Implementation for Backward Adaptive Transform Coding 

An implementation of the invention for backward adaptive transform coding will now be 
described. In adaptive transform coding, adaptation of the transform in the encoder must be 

1 0 accompanied by corresponding adaptation in the decoder. With a typical forward adaptive structure, 
changes to the transform must be encoded and sent along with the encoded data. An alternative is 
to base the transform updates only on data available at the decoder. In this case, the decoder can 
deduce the correct decoding transform by running the same adaptation algorithm as the encoder. 
This is analogous to the use of backward adaptation of prediction filters and quantizer scaling in 

15 adaptive differential pulse-coded modulation (ADPCM), as described in, e.g., R.V. Cox, "Speech 
Coding," The Digital Signal Processing Handbook, Chapter 45, CRC and IEEE Press, pp. 45.1- 
45.19, 1998, which is incorporated by reference herein. 

Of course, quantization is an irreversible reduction in information. For a given accuracy, 
estimation of the source covariance requires more samples when using quantized data. Estimating 

20 moments of the source by moments of the quantized data introduces a bias, which for a Gaussian 
source is described in, e.g., L. Cheded, et al., "The Exact Impact of Amplitude Quantization on 
Multi-Dimensional, High-Order Moments Estimation," Signal Proc, Vol. 39, No. 3, pp. 293-315, 
September 1 994, which is incorporated by reference herein. For uncorrelated components, the bias 
disappears from cross moment estimates, such that the bias from quantization has no effect at 

25 convergence. The primary effect of quantization is to increase the steady-state error when the 
quantization is very coarse. Quantization may speed the initial stages of convergence by amplifying 
small sample values. 

FIG. 5 shows a signal processing device which includes an encoding and adaptation structure 
for backward adaptive transform coding in accordance with the invention. Each of the elements 
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labeled Q represents a quantizer, but the structure is otherwise similar to that of the signal processing 
device of FIG. 2. 

When the transform adaptation is driven exclusively by quantized data that will be 
transmitted to a decoder, it is not necessary to separately encode the transform update. Instead, by 
5 using the same adaptation rules, the decoder can track the transform used in the encoder. 

FIG. 6 shows the results of simulations of a backward adaptive version of Algorithm 3, using 
the same example four-dimensional source X as the previous example. The source in these 
simulations is thus Gaussian with E[xpcj\ = 0.6 1,vl . The initial transform is again the identity 
transform, and the step size is again ji = \i*fy, where \i* is the maximum step size for convergence 

10 in mean. Based on the results depicted in FIG. 4C, the step size ^ = \i*/20 -0.02378 was chosen, 
i.e., y = 20. The quantization is unbounded and uniform with quantization step size A. The curves 
in FIG. 6 correspond to different values of A, the step size of the uniform scalar quantizer used in 
the encoder, and a curve for the no quantization case is included for comparison purposes. Each 
curve represents the average of 50 simulations. 

15 It can be seen from the curves of FIG. 6 that for a reasonable range of quantization step sizes, 

the use of quantized data increases the steady-state error slightly and has little effect on the rate of 
convergence. Several other simulations of this type are presented in the above-cited Ph.D. thesis of 
V.K. Goyal, including cases where the source is time-varying. 

The above-described embodiments of the invention are intended to be illustrative only. 

20 Alternative embodiments of the invention may utilize other signal processing devices, with different 
arrangements and configurations of elements. For example, the invention may utilize reduced- 
parameter forms other than the Givens parameterized form described herein, such as a Householder 
form. The techniques of the invention are also applicable to any desired type of source signal(s). 
Moreover, the invention may be used for a wide variety of different types of signal processing 

25 applications other than those described herein. Furthermore, the above-described processing 
operations of the invention may be implemented in whole or in part using one or more software 
programs. These and numerous other alternative embodiments within the scope of the following 
claims will be apparent to those skilled in the art. 
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